function [pi0,R,T,K] = CalibrateCamera2( Xreal, xImg, PItru_ref )

% [xImg  xmean ymean scalef] = normalizeImgCoord(xImg);
% [Xreal  xmean3 ymean3 zmean3 scalef3] = normalize3DCoord(Xreal);

% This method derives from the cross-product formulation:

% lambda * xi = PI * X   
% xi CROSS lambda xi = 0 = xi CROSS (PI * X ) then writing the two lin.
% indep parts we get:
% 0 = -p2'*X + yi * p3'*X
% 0 = p1'*X - xi*p3'*X

% verify that this is actually zero!!!
vv = cross(xImg, PItru_ref * Xreal );

R = zeros(3,3);
T = zeros(3,1);
K = zeros(3,3);

N = size(xImg,2);
M = zeros(2 * N, 12);
xobs = xImg(1,:);
yobs = xImg(2,:);

% ... this is show stopper for applying this method
% you need to know the depths of the points to calibrate
% but that's what you want to find out!!
% (unless you are willing to measure the lengths physically,
% from 3D point to camera along camera axis)
zref = xImg(3,:);

for n = 1:2:2*N
   k = ceil(n/2);
   M(n,1:4) = zeros(1,4); 
   M(n,5:8) = -zref(k)*Xreal(1:4,k); 
   M(n,9:12) = yobs(k)*Xreal(1:4,k); 
   
   M(n+1,1:4) = zref(k)*Xreal(1:4,k);
   M(n+1,5:8) = zeros(1,4);
   M(n+1,9:12) = -xobs(k)*Xreal(1:4,k);
    
end

[U,S,V] = svd(M); 
pi_star = V(:,12);

pi0 = reshape(pi_star,4,3)';

KR = pi0(1:3,1:3);
KT = pi0(1:3,4);
[q,r] = qr(KR);
Ksq = q*r*(KR)';
K = sqrt(diag(diag(Ksq))); % Assumption: K is diagonal (no skew, offsets are pre-factored)
R = K^-1 * KR;
T = K^-1*KT;
% [K,R] = rq(KR);
% T = K^-1 * KT;

breakhere =1 ;